#ifndef _FIELD_ENERGY_H_
#define _FIELD_ENERGY_H_

#include <map>
#include <numeric>

#include <boost/math/constants/constants.hpp>
#include <CGAL/squared_distance_2.h>
#include <CGAL/centroid.h>
#include <CGAL/Aff_transformation_2.h>

#include "../common/polytope.h"
#include "../common/potential_energy.h"

namespace Field_packing
{
	using namespace Packing_common;

	class Field_energy_system : public System_potential_energy
	{
		typedef Packing_polygon::Point_2 Point_2;
		typedef Packing_polygon::Vector	Vector_2;
		typedef CGAL::Line_2<CGAL::Exact_predicates_inexact_constructions_kernel>	Line_2;
		typedef CGAL::Aff_transformation_2<CGAL::Exact_predicates_inexact_constructions_kernel> Transformation_2;
	public:

		Field_energy_system(const std::vector<Packing_polygon*>& polygon_set, const Container_box& bndbox) ;
		int optimize(double *trans, bool apply_constraints = false) ;
		double system_energy_at(const double *const x);

		// objective function
		double gravity_energy_at(const double *const x) { return 0.0; }
		double molecular_energy_at(const double *const x) ;
		
		~Field_energy_system() { }

	private:
		double point_polygon_squared_distance(const Point_2& p, const Packing_polygon& polygon);
		double point_polygon_squared_distance(const Point_2& p, const Packing_polygon& polygon, unsigned int& edge_idx); 

		void constraint_at(const double *const x, double *constraint_values) ;
		void constraint_grad_at(const double *const x, double *constraint_grad) ;

		void energy_gradient_at(const double *const x, double *gradient) ;
		void molecular_energy_gradient_at(const double *const x, double *gradient) ;

	/*	int callback(const int evalRequestCode, const int n, const int m, const int nnzJ, 
							const int nnzH,	const double * const x,	const double * const lambda, 
							double * const obj, double * const c, double * const objGrad, double * const jac,
							double * const hessian,	double * const hessVector, void * userParams);*/

	private:
		const std::vector<Packing_polygon*>& polygons;
		double rm; // at which the potential reaches minimum
		// parameters used in the attenuation function
		double a;
		double lambda;
	};
}

#endif